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We present a long-wavelength approximation to the Navier-Stokes Cahn-Hilliard 
equations to describe phase separation in thin films. The equations we derive un- 
derscore the coupled behaviour of free-surface variations and phase separation. We 
introduce a repulsive substrate-film interaction potential and analyse the resulting 
fourth-order equations by constructing a Lyapunov functional, which, combined with 
the regularizing repulsive potential, gives rise to a positive lower bound for the free- 
surface height. The value of this lower bound depends on the parameters of the 
problem, a result which we compare with numerical simulations. While the theoret- 
ical lower bound is an obstacle to the rupture of a film that initially is everywhere 
of finite height, it is not sufficiently sharp to represent accurately the parametric 
dependence of the observed dips or 'valleys' in free-surface height. We observe these 
valleys across zones where the concentration of the binary mixture changes sharply, 
indicating the formation of bubbles. Finally, we carry out numerical simulations 
without the repulsive interaction, and find that the film ruptures in finite time, 
while the gradient of the Cahn-Hilliard concentration develops a singularity. 



I. INTRODUCTION 



Below a certain critical temperature, a well-mixed binary fluid spontaneously separates into 
its component parts, forming domains of pure liquid. This process can be characterised by 
the Cahn-Hilliard equation, and numerous studies describe the physics and mathematics of 
phase separation [TJ [2J, |3j HJ E]. In this paper we study phase separation in a thin layer, 
in which the varying free-surface and concentration fields are coupled through a pair of 
nonlinear evolution equations. 

Calm and Hilliard introduced their eponymous equation in [1J to model phase separation 
in a binary alloy. Since then, the model has been used in diverse applications: to describe 
polymeric fluids [6j, fluids with interfacial tension [8], and self-segregating populations 
in biology [9]. An analysis of the Cahn-Hilliard (CH) equation was given by Elliott and 
Zheng [10], where they obtain existence, uniqueness, and regularity results. Several au- 
thors have developed generalisations of the CH equation: a variable-mobility model was 
introduced by Elliott and Garcke [llj, while nonlocal effects were considered by Gajweski 
and Zacharias [12] . These additional features do not qualitatively change the phase sepa- 
ration, and we therefore turn to one mechanism that does: the coupling of a flow field to 
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the Cahn-Hilliard equation [2]. In this case, the Cahn-Hilliard concentration equation is 
modified by an advection term, and the flow field is either prescribed or evolves according 
to some fluid equation. Ding and co-workers [7J provide a derivation of coupled Navier- 
Stokes Cahn-Hilliard (NSCH) equations in which the velocity advects the phase-separating 
concentration field, while concentration gradients modify the velocity through an additional 
stress term in the momentum equation. A similar model has been produced by Lowengrub 
and Truskinowsky [3J. Such models have formed the basis of numerical studies of binary flu- 
ids [13] , while other studies without this feedback term highlight different regimes of phase 
separation under flow [51 HU [15]. Here, the NSCH equations form the starting point for our 
asymptotic analysis. 

As in other applications involving the Navier-Stokes equations, the complexity of the 
problem is reduced when the fluid is spread thinly on a substrate, and the upper vertical 
boundary forms a free surface [TBI ITT] . Then, provided lateral gradients are small compared 
to vertical gradients, a long-wavelength approximation is possible, in which the full equations 
with a moving boundary at the free surface are reduced to a single equation for the free- 
surface height. In the present case, the reduction yields two equations: one for the free 
surface, and one for the Cahn-Hilliard concentration. The resulting thin-film Stokes Cahn- 
Hilliard equations have already been introduced by the authors in [18] . although the focus 
there was on control of phase separation and numerical simulations in three dimensions. Here 
we confine ourselves to the two-dimensional case: we derive the thin-film equations from first 
principles, present analysis of the resulting equations, and highlight the impossibility of film 
rupture, once a regularizing potential is prescribed. We use numerical simulations to show 
that in the absence of this regularizing potential, the film does indeed rupture, an event that 
coincides with the development of a singularity in the concentration gradient. 

Along with the simplification of the problem that thin-film theory provides, there are 
many practical reasons for studying phase separation in thin layers. Thin polymer films 
are used in the fabrication of semiconductor devices, for which detailed knowledge of film 
morphology is required [19] . Other industrial applications of polymer films include paints 
and coatings, which are typically mixtures of polymers. One potential application of the 
thin-film Cahn-Hilliard theory is in self-assembly [2TJ] |2"T| I2"2l |2"3"| 121] . Here molecules (usually 
residing in a thin layer) respond to an energy-minimisation requirement by spontaneously 
forming large-scale structures. Equations of Cahn-Hilliard type have been proposed to 
explain the qualitative features of self-assembly [201 EH]> an d knowledge of variations in the 
film height could enhance these models. Indeed in [18] the authors use the present thin- 
film Cahn-Hilliard model in three dimensions to control phase separation, a useful tool in 
applications where it is necessary for the molecules in the film to form a given structure. 

The mathematical analysis of thin-film equations was given great impetus by Bernis and 
Friedman in [2B]. They focus on the basic thin-film equation, 

dh = _d_ / ^ 
dt dx \ dx 3 J ' 

with no-flux boundary conditions on a line segment, and smooth nonnegative initial con- 
ditions. For n = 1 this equation describes a thin bridge between two masses of fluid in a 
Hele-Shaw cell, for n < 3 it is used in slip models as h — > [27J, while for n = 3 it gives 
the evolution of the free surface of a thin film experiencing capillary forces [IB] . Using a 
decaying free-energy functional, they analyzed Eq. ([T| and obtained results concerning the 
existence, uniqueness, and regularity of solutions, as a function of the exponent n. Only for 
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n > 4 does a classical, smooth solution exist; this fact is established by construction of an 
entropy functional. This paper [26] has inspired other work on the subject [281 12H1 EQ], in 
which the effect of a Van der Waals term on Eq. ([T| is investigated. These works provide 
results concerning regularity, long-time behaviour, and film rupture in the presence of an 
attractive Van der Waals force. More relevant to the present work is the paper by Wieland 
and Garcke [31] , in which a pair of partial differential equations describes the coupled evolu- 
tion of free-surface variations and surfactant concentration. The authors derive the relevant 
equations using the long-wavelength theory, obtain a decaying energy functional, and prove 
results concerning the existence and non-negativity of solutions. 

When the binary fluid forms a thin film on a substrate, we shall show in Sec. [IT] that a 
long-wave approximation simplifies the density- and viscosity-matched Navier-Stokes Cahn- 
Hilliard equations, which reduce to a pair of coupled evolution equations for the free surface 
and concentration. If h(x,t) is the scaled free-surface height, and c(x,t) is the binary fluid 
concentration, then the dimensionless equations take the form 



dt dx ^' dt ^ C dx dx V dx 



(2a) 
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Here C is the capillary number, r measures the strength of coupling between the concentra- 
tion and free-surface variations (backreaction) , and C n is the scaled interfacial thickness - 
sometimes called the Cahn number. Additionally, a is the dimensionless, spatially-varying 
surface tension, and <fi is the body-force potential acting on the film. In this paper we fo- 
cus on a class of potentials that models a repulsive interaction between the film and the 
substrate, thus preventing rupture. This enables us to focus on late-time phase separation, 
which is synonymous with a tendency to equilibrium. However, rupture is in itself an im- 
portant feature in thin-film equations [T6] 128] [29] : we therefore use numerical methods to 
highlight the possibility of rupture in the absence of a repulsive film-substrate interaction. 

The paper is organised as follows. In Sec. [IT] we discuss the Navier-Stokes Cahn-Hilliard 
equation and the scaling laws that facilitate the passage to the long-wavelength equations, 
and we derive Eq. ([2]). In Sec. Ill we perform linear and non- linear analyses of the long- 



wavelength equations. Our non-linear analysis centres on finding a Lyapunov functional for 
a given class of potentials. We derive a priori bounds for a given positive (h > 0) solution, 



and estimate the minimum value of the free-surface height. In Sec. IV we outline a series of 
numerical studies, with and without a regularizing potential, and we discuss the dependence 
of the minimum free-surface height on the problem parameters. Finally, in Sec. |V]we present 
our conclusions. 



II. THE MODEL EQUATIONS: DERIVATION 



In this section we introduce the two-dimensional Navier-Stokes Cahn-Hilliard (NSCH) equa- 
tion set. We focus on the so-called matched case, wherein both components of the binary 
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mixture have the same density and viscosity. We discuss the assumptions underlying the 
long-wavelength approximation. We enumerate the scaling rules necessary to obtain the 
simplified equations. Finally, we arrive at a set of equations that describe phase separation 
in a thin film subject to arbitrary body forces. 

The full NSCH equations describe the coupled effects of phase separation and flow in a 
binary fluid. If the fluids are density- and viscosity matched, then the models described in 
the references E] agree; this is the case we study. If v is the fluid velocity and c is the 
concentration of the mixture, where c = ±1 indicates total segregation, then these fields 
evolve as 

+ u-V« = V- T- -V0, (3a) 
at p 

oc 

— + v ■ Vc = DV 2 (c 3 - c - 7 V 2 c) , (3b) 
V • v = 0, (3c) 

where 

V \dxj dxi) ^ dxi dxj ^ 

is the stress tensor, p is the fluid pressure, is the body potential and p is the constant 
density. The constant v is the kinematic viscosity, v = rj/p, where rj is the dynamic viscosity. 
Additionally, (3 is a constant with units of [Energy] [Mass]" 1 , «/7 is a constant that gives 
the typical width of interdomain transitions, and D a diffusion coefficient with dimensions 
[Length] 2 [Time]" 1 . 

We impose the following boundary conditions (BCs). On the lateral boundaries (re- 
direction), the velocity must satisfy the no-slip condition, while c x and c xxx must vanish too. 
Alternatively, we may enforce periodicity, and demand that the velocity, c, c x , and c xx be 
periodic functions of the lateral coordinate. Finally, if the system has a free surface in the 
vertical or z-direction, then the vertical BCs are the following: 

u = w = c z = c zzz on z = 0, (5a) 

while on the free surface z = h(x,t) they are 

d 

hihjTij = -an, nrfjTij = ——, (5b) 
dh dh . , 

W= m +U ^ (5c) 

hidiC = 0, hidiV 2 c = 0, (5d) 

where ft = (—d x h, 1)/[1 + (d x h) 2 ] 1 / 2 is the unit normal to the surface, t is the unit vector 
tangent to the surface, s is the surface coordinate, a is the surface tension, and k is the 
mean curvature, 

„ d xx h 
k = V • n = 



21 2 



[i + OW] 

these free-surface conditions are standard and are discussed in the review papers [161 EE]- 
This choice of BCs guarantees the conservation of the total mass and volume, 

Mass = / dxdz c(x, z,t), Volume = / dxdz. (6) 

JDom(t) JDom(t) 
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Here Dom(t) represents the time-dependent domain of integration, owing to the variability 
of the free-surface height. Note that in view of the concentration BC (5d), the stress BC (5b) 
and does not contain c(x,t) or its derivatives. 

These equations simplify considerably if the fluid forms a thin layer of mean thickness ho, 
for then the scale of lateral variations I is large compared with the scale of vertical variations 
ho- Specifically, the parameter 5 — ho/£ is small, and after nondimensionalisation of Eq. ^ 
we expand its solution in terms of this parameter, keeping only the lowest-order terms. For 
a review of this method and its applications, see [16] IT?] . For simplicity, we shall work in 
two dimensions, but the generalisation to three dimensions is easily effected [IS] . 

In terms of the small parameter 5, the equations nondimensionalise as follows. The 
diffusion time scale is to = £ 2 /D = h\j {S 2 D) and we choose this to be the unit of time. 
Then the unit of horizontal velocity is uq = £/to = 5D/h so that u = (8D/ho) U , where 
variables in upper case denote dimensionless quantities. Similarly, the vertical velocity is 
w = (5 2 D/ho)W, and the free-surface height is h = h H. The dimensionless coordinates 
are introduced through the equations x = £X, z = hoZ. Finally, for the equations of motion 
to be half-Poiseuille at O (1) (in the absence of the backreaction) we choose p = (r]D/h, 
and = (r/D/hl) $. We also stipulate the following form for the surface tension: 



a = a [1 + <$«/ (X)] , 

where the exponent q is yet to be determined. Using these scaling rules, the dimensionless 
momentum equations are 
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d 2 c d 2 c 
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dX 2 dZ 2 



(9) 



(10) 



V V 

The choice of ordering for the Reynolds number Re allows us to recover half-Poiseuille flow 
at O (1). We delay choosing the ordering of the dimensionless group until we have 

examined the concentration equation, which in nondimensional form is 



dc TT dc „ T dc\ 

WT + U ax + W dz) 

B 2 d 2 



dX 2 



dZ< 



(c 3 - c) - 5 4 C, 



d A c 
n dX~ A 
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Q 8X 2 dZ 2 ' 



(11) 
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where C 2 = 7//1Q. By switching off the backreaction in the momentum equations (cor- 
responding to P'y/Du — > 0), we find the trivial solution to the momentum equations, 
U = W = d x (P + = d z (P + $) = 0, H = 1. The concentration boundary condi- 
tions are then cz = czzz = on Z = 0, 1, which forces cz = so that the Cahn-Hilliard 
equation is simply 

dT dx 2 { } n «9X 4 ' 
To make the lubrication approximation consistent, we take 



5C Q = C n = 6^//h = O (1) 



(12) 



We now carry out a long- wavelength approximation to Eq. (11), writing U = Uq + O (5), 
W = Wq + O (5), c = c + Sci + <5 2 c 2 + • • • . We examine the boundary conditions on c(x, t) 
first. They are n • Vc = n ■ VV 2 c = on Z = 0, H; on Z = these conditions are simply 
dzc = dzzzc = 0, while on Z = H the surface derivatives are determined by the relations 

n-V oc -5 2 H x d x + d z , 

h-VV 2 cc -5 i H x d xxx -5 2 H x d x dzz + S 2 d xx d z + dzzz. 

Thus, the BCs on c are simply dzc = dzzzCo = on Z = 0, H, which forces c = c (X, T). 
Similarly, we find C\ = C\ (X, T), and 



dc 2 
dZ 



d 2 c 2 H x dc 



for any Z G [0, H] 



H dX' dZ 2 H dX 

In the same manner, we derive the results dzzzzc 2 = dzzzz^ = 0. Using these facts, 
Eq. (11) becomes 

dc 



4- U 9C ° 

df + Uo dx 



— (4 - co) - + (3d - 1) ^ - 26?^^ - CS^ 



We now integrate this equation from Z = to H and use the boundary conditions 

<9 3 c 4 
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After rearrangement, the concentration equation becomes 
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is the vertically-averaged velocity. Introducing 



At = c - c 



Cn d_ ( rrdcp 

H dx\ ox 



the thin-film Cahn-Hilliard equation becomes 

dc <9cp_ = 1 d 

dT [ o) dX ~ HdX 



H 



d/i 
dX 



(13) 



We are now able to perform the long- wavelength approximation to Eqs. ([7]) and (|8]). At 
lowest order, Eq. (TSJ) is dz (P + $) = 0, since c = c (X, T), and hence 



P + $ = P 



surf 



surf 



P (X, P(X, T),T) + $ (X, P(X, T), T) . 



We introduce a dimensionless group to measure the strength of the interaction between the 
concentration and velocity fields; we also specify its order of magnitude: 



$ 2 Pl 
Dv 



Oil) 



(14) 



Later on we refer to this quantity as the 'backraction strength', since it is a measure of the 
extent to which concentration gradients feed back into the flow field. Using this dimensionless 
group, Eq. (|7p becomes 



d 2 U 
dZ 2 



^(P surf + $ surf )+r— I — 



_l_ r 



dc Q d 2 c 2 
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Using dzzC2 = (H x /H) (dc /dX) this becomes 
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HdX 



H 



dc 
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(15) 



At lowest order, the BC (5b) reduces to 

dU 



<9£ 



dZ dX 



on Z = H, 



(16) 



which combined with Eq. (15) yields the relation 

dU 



dZ 



<9£ (8 

dx + [z ~ H) \dx (Psurf + $surf) 



r d 
HdX 



H 



dc 
dX 



Here S is the dimensionless, spatially-varying component of the surface tension. Making use 
of the BC Uq = on Z = and integrating again, we obtain the result 



U (X, Z, T) = Z— + {\Z 2 - HZ) { — (P surf + $ surf ) ' 



\dX 



HdX 



H 
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The vertically-averaged velocity is therefore 

<9£ 1 A d f ld 2 H 



(Uo) 



* H 'dX 



1 Tj2 

* H \dX 



r d 

cdX 2+ swl + HdX 
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(18) 



where we used the standard Laplace- Young free-surface boundary condition to eliminate 
the pressure, and 

v P D „, t , (ig) 



C 



or 



h a 5 2 

Finally, by integrating the continuity equation in the Z-direction, we obtain, in a standard 
manner, an equation for free-surface variations, 



OH d 
df + d~X 



(H(U o ))=0. 



(20) 



The inclusion of the surface-tension terms requires further elucidation. In the long-wave 



limit, the normal-stress condition reduces to p\ 



-ahrr, which in dimensionless form is 



P 



' vpD 



[l + 6«f(x)]H xx . 



By promoting the constant C = vpDj (h <T 5 2 ) to O (1), and by taking q > 1, this equation 



reduces to P = —C X H XX , as in Eq. (18). Similarly, the normal-stress condition reduces to 



vp(du/dz) h = da/dx, which in dimensionless form reads 

dU _ a h 5 q df 
dZ ~ vpD dX' 

Taking q = 2 gives a contribution to the shear-stress balance at lowest order, dU/dZ 



dYi/dX, S = f (X) /C, as in Eq. (16). Going to higher exponents q > 2 suppresses this 
contribution. 

Let us assemble our results, restoring the lower-case fonts and omitting ornamentation 

(21a) 



over the constants. The height equation (|20|) becomes 
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while the concentration equation (13) becomes 
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and where we have the nondimensional constants 

The boundary conditions are inherited from the full Navier-Stokes Cahn-Hilliard equations: 
Since J contains a depth-averaged velocity, we write it as J := hu. Thus, the concentration c, 
the chemical potential /i, and the flux J are either periodic functions in the lateral direction, 



or satisfy the no-flux conditions c x = [i x — J = on the lateral boundaries. Equations (21 ), 
together with the boundary conditions described, are the thin-film NSCH equations. The 
integral quantities defined in Eq. (|6| are manifestly conserved, while the free surface and 
concentration are coupled. The term (h 2 /2)a x in the flux J represents a driving force, 
which can be externally prescribed, or a function of the concentration c. In either case, the 
inclusion of this term can have a substantial effect on the behaviour of the system. For the 
rest of this study, this term is set to zero; its inclusion is discussed elsewhere by the present 
authors [TBI , and by others [32J. 



In view of the severe constraint C n = 5^ /ho = 0(1), some discussion about the ap- 



plicability of Eqs. (21) to real systems is warranted. This constraint is the condition that 
the mean thickness of the film be much smaller than the transition-layer thickness. In ex- 
periments involving the smallest film thicknesses attainable (10~ 8 m) (33], this condition is 
naturally satisfied. Furthermore, in certain situations far from this limiting case, variations 
in the domain structure in the vertical direction are suppressed, and a system of equations 



with no vertical (z-) dependence, such as Eqs. (21), is appropriate. This kind of situation 
arises when external effects such as the air-fluid and fluid-substrate interactions do not pre- 
fer one binary fluid component or another; hence, the dimensionality of the film is reduced, 



and the balance laws implied by Eqs. (21) are applicable. 



III. THE MODEL EQUATIONS: ANALYSIS 



The choice of potential (f>o determines the behaviour of solutions. In this section, we 



perform a linear-stability analysis on the model equations (21) and identify the pattern 



formation mechanism. We also develop results for the non-linear regime using the theory of 
a priori bounds. These are bounds on norms of the solution (h, c) that are obtained without 
assuming any prior knowledge of the solution. Throughout this section, the driving force 
resulting from surface-tension gradients is set to zero. 

The first step in our analytical study is to find the circumstances under which the con- 
stant state (h , Cq) is unstable to a small-amplitude, initial perturbation (Sh ,Sco). This 
perturbation evolves in time to a state (d~h,d~c) (x,t), which satisfies the linearized version 
of equations (21). By writing down a wave ansatz (5h,5c) oc e lkx , we obtain an eigenvalue 
equation, 



d_(6h 
dt V 5c 

with eigenvalues 







-k 2 (3c 2 - 1) - Clk A 



5h 
8c 



A, 



hlk 2 



C 



+ 0o (ho 



(3c 2 - 1) k 2 - C 2 n k A 



(23) 
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Thus, there are two routes to instability. The system can become unstable as a result of 
substrate-film interactions if <p' (ho) < 0. Such an interaction will often lead to rupture |16j . 
If this route is suppressed, then film rupture may be prevented, but the second route to 
instability is also relevant. This is accessible when cq is in the spinodal range | cq | < l/\/3- 
Thus, even when the first route to instability is not accessible, a critical mixed state will 
phase separate in a manner similar to the classical Cahn-Hilliard fluid, as described in Sec.|T} 
While the linear analysis is helpful to describe early-stage evolution, it sheds no light on 



the behaviour at later times. We therefore turn to the non-linear analysis of the problem ( 21 ) . 
The non-linear analysis centres on finding bounds for a given solution (h,c). To do this, 
it is necessary to construct a Lyapunov functional, that is, a non-negative, non-increasing 
functional of the solution (h,c). It is certainly the case that we can find a non- increasing 
functional based on the solution pair (h, c), which is a kind of energy for the problem: 

Proposition 1 (Existence of a decreasing functional) Given a smooth solution (h, c) 
to the equations (21), positive in the sense that h(x,t) > 0, and a continuous potential 
function 4>o, then the functional 



T [h, c] 



dx 



J_ fdh\ 
2C\dx) 



[s] ds 



r<2 



dx h 



!( C 2 -i) 2 + 



dx J 



is non-increasing, T < 0. 



(24) 



The proof of this claim is readily obtained by a straightforward time-differentiation of T [h, c] 
and application of the equations (21), together with the no-flux or periodic boundary con- 
ditions. We find, 

-j\^ + hrf)<0. (25) 



We build on this result by focussing on the following class of potential whose anti-derivative 
is positive: 



>i « 



V) ds' > 0, 



< s < a, 



(26) 



where a is an arbitrary reference height. Using Prop. [I] and the condition (26), we obtain a 
Lyapunov functional for the positive solution (h,c), h > 0: 

Proposition 2 (Existence of a positive Lyapunov functional) For a smooth solution 



(h, c) to the equations (21 ) ; positive in the sense that h (x, t) > 0, and for a potential function 



</>o with positive anti- derivative, there is an associated Lyapunov functional. 

To verify this claim, it suffices to note that since <f>\ > for the class of potential functions 
under consideration, all terms in the functional T [h, c] (Eq. (J24l)) are positive, and thus T 
is a positive, non-increasing function of time, i.e. a Lyapunov functional. 

The boundedness result T it) < T (0) provides a regularity condition on the height 
h (x,t), although this is valid only in a single spatial dimension. 

Proposition 3 (Holder continuity of h(x, •)) If (h,c) is a smooth, positive solution to 
the equations ( pZlj ), in the sense that h(x,t) > 0, and if the potential function 0o has a 
positive anti-derivative, then h(x,-) is Holder continuous, with time-independent Holder 
constant kn- 
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Proposition [3] follows from the a priori bounds 

^£d x h 2 x <F(t)<F(0), 
and from the use of Holder's inequality on the following string of relations: 



(27) 



\h (x 2 ) - h (xi) 



/•X2 PX2 I rL 

I h x dx < I dx\h x \ < \x 2 — Xi\ 1 / 2 \l I dx h? x < ku\x 2 — ^l] 1 ^ 2 , 



where k# = (0) is the time-independent Holder constant. As an immediate corollary 

of this result, we obtain an upper bound on the height field: 

Proposition 4 (An upper bound on the height field) If (h, c) is a smooth, positive 
solution to the equations (21), in the sense that h(x,t) > 0, and if the potential function 4>o 
has a positive anti-derivative, then h(x,-) is bounded above. 

Since the free energy contains a term in the L 2 -norm of h l l 2 c x , a similar result exists for 
the concentration field, provided h > everywhere. Loss of control over the minimum value 
of h therefore implies loss of control over the concentration gradient. This suggests that 
blowup of gradients and film rupture are related, a claim which we demonstrate numerically 
in Sec. |IVB Such extreme events are avoided when a repulsive film-substrate interaction is 
present, in which case a positive lower bound on h exists; it is to that result that we now 
turn. 

Given the form of Eqs. (21), regularity of a given solution (h, c) is guaranteed only when 
a lower bound on h is obtained, in addition to the upper bounds just provided. To derive 
such a result, we first specialize to the potential 



A G 
0O = -2^ 

from which a more general result will follow. 



G > 0. 



(28) 



Proposition 5 (No-rupture condition for the potential in Eq. (28)) If (h,c) is a 
smooth, positive solution to the equations (21), in the sense that h(x,t) > 0, and if the po- 
tential function 4>o has the form given by Eq. (28), then there is an a priori, time-independent 
lower bound on h. 



Note first of all that the potential (28) has a positive anti-derivative, <f>\ (s) = Gs~ 2 , where 
the reference height a in Eq. (26) is set to a = oo. Thus, there is a Lyapunov functional for 
the solution (h,c), and hence, 



j dxGh- 2 = dxfaih) <F{t) <T{$). 
Jo Jo 



Using the Holder continuity of h, 

h 2 < h 2 min + 2h min k H L^ 2 + k 2 H \x- x min \ < h 2 min + 2h min k H L 1 / 2 + k 2 H x, 
where ku = ^2CT (0) is the Holder constant for h. Using results (29) and (30), 

dx 



G 



h 2 mbl + 2h min k H LV 2 + k 2 H x 



<^(0). 



(29) 



(30) 



(31) 
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By integrating this equation, we arrive at the relation 



, / k 2 L \ .F(0U 2 



h 2 nin + 2k H L^h min J ~ G 

Let us examine the properties of the function log [l + k\Lj (s 2 + 2fcj/L 1 ' 2 ,s)] . It tends 
to infinity as s — > 0, and tends to zero as s — ► oo. It is also monotone-decreasing over 
s e (0, oo). Thus, the equation log [l + k 2 H L/ {s 2 + 2k H L l l 2 s)\ = J 7 {0) k 2 H /G has precisely 



one positive root for J 7 (0) ^ 0, which we call s*. To satisfy the inequality (J32J) , it must be 
the case that 

^min ^ *5* ^ 0- 



For the potential 0o = — (G/2) s (Eq. (28)), the root s* can be obtained explicitly: 

"mm d. — n-H-^ 



exp (JF(0) fc^/G) 



1 



> 0. 



This completes the proof of the no- rupture condition for the potential (28) 



To arrive at a no-rupture condition for a general potential, we introduce the function 

2 (s):= I dx 0i (s + k H x 1/2 ) . (33) 
Jo 

Based on Prop. |5j we write down sufficient conditions on 0i and 02 for the existence of a 
positive lower bound on h: 

Proposition 6 (A sufficient condition to avoid rupture) // (h, c) is a smooth, posi- 



tive solution to the equations (21), in the sense that h(x,t) > 0, if the anti- derivative of the 
potential function 0o is a positive, non-increasing function, and moreover, if 02 satisfies the 
conditions 

lim0 2 (s) = oo, 

s^O 

lim 02 (s) < 0, (34) 

s— >oo 

then a positive a priori lower bound on h (x,t) exists, independent of time. 

The proof of Prop. [6] is in the same spirit as that of Prop. [5j Given the positivity of the 
anti-derivative 0i, a Lyapunov exponent exists (Prop. [TJ, and thus we have the bound 
J dx<p\ (h) < T (0). Using the Holder continuity of h (Prop. [3]), and the condition that the 
0i should be a non-increasing function, 

h(x,t) < h min + k H x 1/2 , 
& (MM)) > 0i (^min + k H x 1/2 ) . 

Using the bound f dxcpi {h) < J 7 (0), 

02 (/imin) = / dx 01 (/l min + k H X 1/2 ) < T (0) . 
JO 
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Given the conditions (34), there exists at least one solution to the equation 2 (s) = F (0), 
for J 7 (0) > 0. Indeed, since 0! is non-increasing, so too is 02, and thus this equation has 
precisely one solution, which we call s*. Then, for the condition 2 (^min) > J 7 (0) to be 
satisfied on the interval (0, oo), we must take 



^min ^> 0. 



Finally, using this theory, we investigate the potential 

G 



ns 



n+l ' 



(35) 



and calculate the n-values for which a no-rupture condition can be found. 



Proposition 7 (Conditions on the potential (35) to avoid rupture) // (h, c) is a 

smooth, positive solution to the equations (21), in the sense that h(x,t) > 0, and if the 
potential function O is given by Eq. (35), then a no-rupture condition is guaranteed to hold 
for n > 2. 



The proof of Prop. |7| follows by a straightforward evaluation of the integral (33). The case 
n = 2 is covered by Prop. [5j Thus, we focus on the case n ^ 2, where the integral 02 (s) has 
the value 



2G 



(n - 2) (n - 1) k% a 



n-2 



a 



a 



n-l 



n 



1 



a 



a 



s/k H , (36) 



and where we have set L — 1. For n > 2 the conditions (34) hold: lim^o 02 (s) = oo, 
lim^oo 02 (s) = 0; note also that 2 (s) is non-increasing. Thus, the equation 2 (s) = T (0) 
has exactly one positive root s*, and this serves as a lower bound on h, h min > > 0. 



Note that the relation given in Eq. (36) fails to satisfy conditions (34) when n < 2. Thus, 



+ 

oa n 

■a 
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FIG. 1: A plot of the integral Jq dx (s + x 1 / 2 ) ™ as a function of s. For n > 2 the integral diverges 
as s — > and tends to zero as s — > oo, allowing for the existence of a positive root for the equation 



f Q dx (s + x 1 ' 2 ) = const. > 0. 



a sufficient condition for the solution (h, c) to possess a no-rupture condition is for the 
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potential 0o to have the form given in Eq. (35), with n > 2. This analysis is also described 
schematically in Fig. [IJ where a plot of the integral J Q L dx (s + x 1 ^ 2 ) n is shown as a function 
of s. For n > 2 the integral diverges as s — > and tends to zero as s — > oo, which provides 
for the existence of a positive root for the equation J Q L dx (s + x 1 / 2 ) ™ = const. > 0. 

Two questions arise from our results. The first involves the form of the potential used in 
the derivation of a no-rupture condition: can we replace the repulsive power-law form with a 
more general function and still obtain a no-rupture condition? The answer to this question 
comes readily through the observation that our construction of a positive lower bound for h 
relies only on the surface-tension and body-force components of the free energy, namely 

f L r i 

F\ = ^ dx —h 2 x + 0i (h) 

Thus, Propositions |5]-[7] can be viewed in the context of the PDE theory for a single variable 
(the free-surface height), where a no-rupture condition exists [31], both for the power-law 
type potential considered here and for the Lennard- Jones potential 0o — —Gs~ ni + Bs~ n2 , 
where G and B are positive constants and n\ > max (712, min (1 + 2n2, 3)) [31]. Hence, our 
regularity results are generalizable to this wider class of potential. 

Having weakened the sufficient condition for the avoidance of rupture, it is reasonable 
to ask, is the presence of a suitable potential even a necessary condition? This question 
is motivated by the single- variable theory for the free-surface height, where under certain 
conditions an entropy functional facilitates the construction of an a priori lower bound 
on the height h(x,t) [2B]. The entropy is obtained through the following steps, which we 
describe for the single- variable case 3h t + {h n [(h xxx /C) — (f) x ]} x = 0: 

1. Identify the power of h that is a factor in the flux J; this is the mobility, m . For the 
single- variable case, m (h) = h n . 

2. Obtain the function m\ (s) = f s ds' f s ds" [m (s")] _1 . 

3. The entropy is then defined as S = dxrrii (h); for the single- variable case, this is 
S = dx h~ n+2 (we have omitted the unimportant constant of proportionality). 

For the single- variable case, 

S(t) + — / dt' / dxh 2 xx + \ I dt' I dx h 2 x (t)' (h) = S (0) > 0, (37) 
3^ Jo Jo Jo Jo 



and the entropy is a non-increasing functional of t, dS/dt < 0. Setting O = in Eq. (37) 
gives the relation 

f L dx 2 f* f L 

/ FI 1 ,m»-2 + / dt l dxh 2 xx = S(0)>0. (38) 
Jo [h(x,t)\ 3G J J 



For n > 4, Holder continuity combined with the bound in Eq. (38) enables the construction 
of a pointwise lower bound on h. When n = 3 (the case considered in this paper), no such 
pointwise bound exists; then, estimates for the entropy based on the Holder continuity of 
h are non-singular in the minimum height. However, a more definitive obstacle than this 



exists when we seek to construct an entropy functional for the system (21), namely, that the 



entropy functional implied by Eqs. (21 ) fails to be a non-increasing function of time: 
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Proposition 8 (The time-derivative of the entropy associated with Eqs. (21) is 



not sign definite) Given a smooth, positive solution (h,c) to the equations (21), in the 
sense that h(x,t) > 0, the rate of change of the functional J Q L dxh a , a < 0, is not sign- 
definite. 

This result is established by direct computation: 



4 / dxh a =~ a{a „ i} r dx (h a+i h x ) h 

dtJo 3C J K x >* 



" ( " 1} fL dx h a+x h 2 J G (h) - Ta (< * 1} [ L dx (h a h x ) x hcl 



3 J XTUK ' 3 

No value of a can give the integral in the last term, namely 



j dx (h a h x ) x hc 2 x = dx [ah a c 2 x h 2 x + h a+1 h xx c 2 x ] 
Jo Jo 

a definite sign. Hence, the entropy functional fails to be non-increasing, and therefore cannot 
be used to construct a lower bound on h. The existence of a suitable potential function in 
the evolution equation for h is thus a necessary and sufficient condition for the avoidance of 



rupture. Indeed, the existence of the backreaction in the equations (21 ), and its consequences 
for the entropy functional, suggest that it promotes rupture while the regularizing potential 
inhibits it. To examine this effect further, we turn to numerical simulations. 



IV. NUMERICAL RESULTS 



In this section we obtain numerical solutions of the equations (21), with and without the 
regularizing potential. Our numerical method is described and validated in Appendix A. 
Where present, the regularizing potential is assigned the form <ft = —G/h 3 , G > 0, which 



satisfies the no-rupture condition given in Sec. Ill In the case where the no-rupture condition 
is satisfied, we study solutions that tend towards an equilibrium. We then characterize 
this equilibrium through the solution of a boundary-value problem. The tendency towards 



equilibrium is in agreement with the predictions of Sec. Ill, where a linearly unstable state 



evolves to reduce the energy functional (24), such that domains of concentration form, 



separated by transition zones, across which the height of the film decreases markedly, forming 
'valleys'. For the case where the no- rupture condition is not automatically satisfied (G = 0), 
we perform numerical studies that suggest that rupture does indeed occur, and coincides 
with the development of a finite-time singularity. Throughout this section, the driving force 
resulting from surface-tension gradients is set to zero. 



A. Numerical studies with a regularizing Van der Waals potential 



We perform numerical simulations of the full equations (21 ), with initial data comprising 
a perturbation away from the unstable steady state (h, c) = (1, 0): 



h (x, 0) = 1, c (x, 0) = 0.01 sin [5 (2tt/L) x] 
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We fix the parameters G = C = 1, and C n = 0.1. The C n - value is small enough such 
that the transition width of domains is much smaller than the size of the computational 
domain L = 2tt. We vary the parameter r between 0.1 and 1 to examine the effects of 
the backreaction strength. The calculations are carried out on a periodic domain, with 
iV = 256 gridpoints and a timestep At = 10~ 4 . The early-stage growth in the concentration 
is governed by the linear theory of Eq. (23). The free energy is a decreasing function of time 
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FIG. 2: Numerical simulation of the decay of the free energy over time, (a) r = 0.1; (b) r = 1 (the 
other parameters are kept constant: C = G = 1, C n = 0.1). There is a very small, local increase 
in the free energy in case (b) due to numerical error. This is small however, and is in contrast to 
the otherwise decreasing trend in the free energy. The energy decreases particularly sharply when 
domains merge: the sharp drops correspond precisely to the domain-merging events in Fig. [3] 

(Fig. m. The average height is conserved exactly by the numerical simulation, while the 
mass Jchdx deviates from its initial value of zero within a range ±10~ 4 for the r = 0.1 case, 
and ±10 -2 for the r = 1 case. In the course of this evolution, the amplitude of the initial 
sinusoidal concentration field grows transiently in time; later on the positive and negative 
'domains' formed by each half-period of the sine function merge to form larger domains, 
each of larger amplitude. At the same time, the free-surface height decreases in value at the 
borders of these domains, forming valleys. Eventually, and as a consequence of the energy- 



minimization principle (25), only a single domain remains. This evolution is best described 



visually, as in Fig. |3j where spacetime plots of c and h show eventual coalescence into a 
pair of opposite-signed c-domains. The domain coalescence happens more rapidly for the 
r = 1 case, compared to the r = 0.1 case. This shows that the coupling of the free-surface 
height to the Cahn-Hilliard concentration far from arresting the domain coarsening, actually 
enhances it. 

The free surface and concentration evolve to an equilibrium state where the salient feature 
is the formation of domains (intervals where c ~ ±1) that are separated by smooth transition 
zones, across which the free surface dips below its average value. We therefore shift focus 
to this state, obtained by setting \i = constant, u = in Eq. (21): 



1 d 2 h 
Cdx^ 



CiG 



1 



1 



d*c 
dx 2 



1 dh dc 

h dx dx ' 



0c 
dx 



(39a) 



(39b) 
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FIG. 3: Temporal evolution of the free-surface and concentration fields. Across the top the backre- 
action strength is set to r = 0.1. Subfigure (a) shows the concentration for this case; (b) shows the 
free-surface height. Across the bottom the backreaction strength has been increased to r = 1.0. 
Subfigure (c) shows the corresponding concentration; (d) shows the free-surface height. The other 
parameters are kept constant: C = G = 1, C n = 0.1. The domains coalesce until only a pair of 
opposite-signed domains remain. The coalescence is more rapid for the r = 1 case, compared to 
the r = 0.1 case. 



where we have enforced the boundary conditions h (±oo) = 1, /i (±oo) = and have rescaled 
lengths by C n . For the case C = oo and p = rjC\G <C 1, Eqs. (39) have an asymptotic 
solution. We find the /i-equation 



1 + P 



^-) 2 + i(|) : 



-1/3 



Hence, the c-equation is 



^_i + lp(c 2 -i) 2 + M!y 



9x2 i+1p(c 2 -i) 2 +MD 2 

For small p, the solution is c = tanh (x/\/2) + O (p) and hence 



(c 3 -c). 



h = 1 - |psech 4 f-^j + O (p 2 ) , p < 1. 



(40) 



(41) 



(42) 
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Thus, in this limiting case, the height profile is approximately constant (h = 1) except in 
the transition region of the concentration field, where it dips. 

The results for the case with finite surface tension are qualitatively similar. Here, two 



parameters characterize the problem: since Eq. (39a) can be multiplied across by C, there 
are precisely two dimensionless groups, CC\G and rC . In Fig. [4] we present numerical 
solutions exhibiting the dependence of the solutions on these parameters. As before, the 
height field possesses peaks and valleys, where the valleys occur in the transition region of 
concentration. While the valley increases in depth for large r or small G, rupture never 



takes place, as guaranteed by the analysis of Sec. Ill The repulsive Van der Waals potential 
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FIG. 4: Solutions of the thin-film equations obtained by solving the boundary-value problem (39). 
Across the top: the effect on the equilibrium solutions of varying the backreaction strength, for 
parameter values C = C^G = 1 and r = 0.1,1,10,50. In (a) the valley deepens with increasing 
r although the film never ruptures, while in (b) the front steepens with increasing r. Figs, (a) 
and (b) are taken from 6 Naraigh and Thiffeault [18J. Across the bottom: the effect on the 
equilibrium solutions of varying the strength of the regularizing potential. The parameter values 
are C = C„ = r = 1 and G = 0.001,0.01,0.1,1. In (c) the valley deepens with decreasing G 
although the film never ruptures, while in (d) the front steepens with decreasing G. 



therefore has a regularizing effect on the solutions. Indeed, the formation of the valley in 
the height field has the physical interpretation of a balance between the Van der Waals 
and backreaction effects. From Fig. [5] we see that the backreaction force, which, through 
Eq. (21) we identify as F cap = —rh~ 1 d x [ h (d x c) 2 ~\ , is of opposite sign to the Van der Waals 
force -Fydw = Gd x h~ 3 . Now in Sec. |III| we showed how the Van der Waals force inhibits 
rupture; hence, _F cap must promote it. The depth of the valley in the height field is therefore 
selected through a balance between rupture-preventing and rupture-promoting effects. 
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FIG. 5: Boundary- value solution. A plot of the forces F cap ana -r v dw 



and -F v d\v for C = = G = 1 and 
r = 50: they have opposite sign; hence, the regularizing potential opposes the rupture-inducing 
tendency of the backreaction. 



Finally, we compare the results for the minimum free-surface height implied by the so- 
lution of the boundary value problem (BVP) with the theoretical lower bound obtained 
in Sec. 



Sec. [Ill 



III In terms of the physical parameters of the system, the no-rupture condition of 



is 



^2CL(F + FiG) 



e iCG- 1 (F +FxGf 



1 >0, 



where F\ = \ f n dx [h (x, 0)] 2 ^ 0, and F = JF(0) — F\. The function s* (G,C) has no 
explicit r-dependence: although F depends on r, it is possible to find initial data to remove 
this dependence. We show a representative plot of s* (G, C) in Fig. [6] (a), while in (b) we show 
a plot of /i m in as a function of G, obtained from the solution of the BVP. Now although these 
two figures represent solutions to the model equations for different boundary conditions, a 
comparison between them is warranted, especially at a domain boundary, where the film 
thinning is induced by entirely local effects. The shape of the two bounds in Figs. [6] (a) 
and (b) is different. Since the bound in Fig. [6] (b) is obtained from numerical simulations, 
and is intuitively correct, we conclude that it has the correct shape and that the bound of 
Fig. |6j while mathematically indispensable, is not sharp enough to be useful in determining 
the parametric dependence of the dip in free-surface height. 



B. Numerical studies without a regularizing Van der Waals potential: a study of 

rupture 



We perform numerical simulations of the full equations (21), with the following finite- 
amplitude initial data: 

h(x,0) = 1 + 0.1 sin [3 (2tt/L) x + tt] , c (x, 0) = 0.5 sin [3 (2tt/L) x) . 

We use the parameters r = 2, G — 0, C — 1, C n = 0.1, and L = 2ir. The calculations are 
carried out on a periodic domain, with 512, and 1024 gridpoints and a timestep At = 10 -6 . 
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(a) (b) 



FIG. 6: (a) A typical plot of s* (G, C) for Fq = F\ = \ and C = 1. This theoretical lower bound has 
a different shape from that in (b), which is obtained from a solution of the BVP, with C = r = 1. 
This suggests that while s* (G, C) plays an important role in the analysis of the model equations, 
it does not capture the physics of film thinning. 



The regularizing Van der Waals force is no longer present, and thus the estimates of Sec. Ill 
no longer apply. We therefore examine the possibility that the film will rupture in finite 
time. For the parameter values chosen, rupture does indeed occur in finite time, as evidenced 
by Fig. [7} In the numerical simulation, rupture is only hastened by grid refinement: this 
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FIG. 7: (a) Temporal evolution of the minimum free-surface height. The rupture is hastened by 
grid refinement; (b) our numerical simulation also captures the decay of the free energy; (c) the 
rupture coincides with a finite-time singularity, wherein the derivative c x diverges at the point 
where h touches down. 



indicates that the effect does not disappear with an increase in resolution, but is nevertheless 
difficult to capture precisely. The simulation also decreases the free energy. This is consistent 



with the free-energy decay law derived in Sec. |III[ which relies only on the fact that h 
should be positive, and holds even for zero Van der Waals forces. Thus, we are satisfied 
that the rupture is accurately described by the numerical simulation, and is not simply an 
artefact. The evolution towards rupture is shown again in Fig. [8] The figure demonstrates 
yet again that rupture is induced in the transition zones of concentration. Our understanding 
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(a) (b) 

FIG. 8: No regularizing potential and large backreaction effect, r = 2. Temporal evolution of (a) 
the concentration; (b) the free-surface height. The height touches down to zero in finite time. A 
singularity develops in the equations and the gradient of the concentration diverges in the transition 
zone. 



of rupture is strengthened by a further examination of the equilibrium case described by 
Eqs. (39). For G = 0, the /i-equation reads 



1 d 2 h 



(43) 



As a boundary- value problem, with boundary conditions h = 1 as x — > ±oo and c = ±1, 
c x = as x — > ±oo, Eq. (43) has no solution, since h xx > everywhere is not compatible 



with a bounded solution. In the language of dynamical systems, the proposed boundary 
conditions for the solution pair (h,h x ) imply a homoclinic orbit, which is impossible if h xx 
is positive everywhere. Thus, the development of a finite-time singularity is consistent with 
the non-existence of a time-independent solution of the (h, c) equation pair. 

Our numerical study has demonstrated the sharp difference in the two cases wherein the 
repulsive Van der Waals force is either included or neglected. In the case where this force is 
neglected, the film ruptures in finite time, an event that is accompanied by the development 
of a singularity in the derivative in the concentration. While interesting mathematically, 
this is undesirable from a physical point of view. Consistent with the analysis of Sec. |III[ 
the inclusion of the Van der Waals term prevents rupture, and enables the development 
of an equilibrium state. Note finally that simulations involving two lateral directions have 
elsewhere been carried out by the authors [18], and the qualitative features are similar to 
those obtained here. 



V. CONCLUSIONS 

Starting from the Navier-Stokes Cahn-Hilliard equations, we have derived a pair of nonlin- 
ear parabolic PDEs that describe the coupled effects of phase separation and free-surface 
variations in a thin film of binary liquid. Since we are interested in the long-time outcome 
of the phase separation, we have focussed on liquids that experience a repulsive Van der 
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Waals force, which tends to inhibit film rupture. Using physical intuition, we identified a 
decaying energy functional that facilitated analysis of the equations. Based on this decaying 
energy functional, we have developed a series of a priori estimates for positive solutions 
(h,c), h > to the model equations (21). The Holder continuity of h obtained through the 
decaying energy gives rise to a positive lower bound on the height h, valid for a general fam- 
ily of repulsive potentials. These estimates are valid not only in the a priori sense described 
here, but also as a means of demonstrating the existence of regular solutions to Eqs. (21), 
given appropriate initial data [26| 135]. 



We carried out one-dimensional numerical simulations of the equations (21) and found 
that the free-surface height and concentration tend to an equilibrium state. The concentra- 
tion forms domains; that is, extended regions where c ~ ±1. The domains are separated 
by narrow zones where the concentration smoothly transitions between the limiting values 
±1. At the transition zones, the free surface dips below its mean value to form a 'valley', 
a feature of binary thin-film behaviour that is observed in experiments. To study the val- 
ley depth as a function of the problem parameters, we focussed on solving the equilibrium 
version of Eq. (21) as a boundary- value problem. This simplification is carried out with- 
out much loss of generality, since our numerical simulations indicate that the system tends 
asymptotically to such a state. We have shown that the valley becomes shallower upon 
increasing the strength of the repulsive Van der Waals force, while it deepens when the 
backreaction strength is increased. The film-thinning tendency of the backreaction has been 
observed experimentally (SHI EH EE] • In the limit of zero repulsive Van der Waals forces, the 
solution of the boundary-value problem implies that the film ruptures, and our temporal 
numerical simulations confirm this. Indeed, we have demonstrated that the film ruptures 
in finite time; simultaneously, the derivative of the concentration becomes singular. Since 
such singularities are undesirable from a physical point of view, this result underscores the 
importance of including Van der Waals forces in studies like this one. 
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APPENDIX A 



In this section we outline a numerical technique for the equations 
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dx J 



CI d 



c — c — 



h dx V dx 



h 



dc 



(Alb) 



(Ale) 



We impose periodic boundary conditions on the solution (h, c) and its derivatives. The 
solution (h, c) is discretized on a regular spatial grid such that the vector pair (h (t) ,c(t)) G 
represents the discretized solution at time t. The derivatives are approximated as 



DiV 



X 



centred finite differences with the periodic boundary conditions taken into account. Thus, 
the derivatives are reduced to matrix operators Vj on R N (here the subscript 'j' denotes 
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the order of the derivative) . The solution is marched forwards in time using a semi- implicit 
Euler algorithm, 



h n+1 - h Tl 
At 

c n+l _ c r< 
At 



1 

3C 



v, [(h n y 3 ■ {v 3 h n+1 )} + sz, 



-Cl [V 4 c n+1 ] + S%, 



(A2) 



where the 'dot' is used here to denote pointwise multiplication, x ■ y = (xxyi, ■ ■ ■ ,xn]Jn), 



x a = (x a 



5 X N)l 



and where the 'source' terms Sh and S c are defined as follows: 



S r . = 



+ r (ft" 1 ) ■ [V, (h-(V lC y 2 )] 



u 



u ■ {V x c) + (h- 1 ) ■ {Vxh) ■ (X>i/x) + V 2 [c 3 -c-Cl (h~ l ) ■ (V x h) ■ (Pic 
1 V 3 h + V l( f> + r (h- 1 ) ■ [T>i (h ■ (Pic)' 2 )] 



C 



c 3 - c - Cl (h- 1 ) ■ {V x h) • (Pic) - ClV 2 c. 



Equations (A2) are re- written as 



h n+l = h n + AtS n^ 

[1 + AtClVi] c n+1 = c n + At 5 6 n , 



(A3) 

(A4) 
(A5) 



a linear problem that is solvable for (h n+1 ,c n+1 ) by matrix inversion. The fa-equation 



manifestly conserves the sum Yli=i since Yli=i Sj=i ^^M v i = 0' f° r an y vec t° r v £ 
Other semi-implicit numerical schemes that only approximate this conservation law can fail 

1.4r 




FIG. 9: Comparison with the work of Burelbach et al. for the rupture of a single-component fluid 
under the influence of an attractive Van der Waals potential. The rupture happens in finite time 
and is calculated here as tn = 4.093, for simulation parameters At = 10~ 5 and N = 100. 



near rupture (h — > 0). One such technique is the otherwise successful method of Kondic [39J, 
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which replaces the term (h n )' 3 -(V 3 h n+1 ) 
iteration at each timestep. 



with (h 



n+l\ 



■(V 3 h 



and thus involves a Newton 



Now although the implied conservation law dt f Q chdx is not manifest in the c- 



equation (A5), we have verified that a sufficiently small stepsize and gridsize guarantees 



its conservation in practice. The implicit step in Eq. (A5) is particularly fast because the 
matrix 1 + AtC^D^ need only be inverted once. This implicit treatment of the high-order 



derivatives in Eqs. (A4)-(A5) also ensures numerical stability for large timesteps that would 



otherwise cause numerical blowup. 

We verify the correctness of our numerical scheme by comparing it against a set of well- 
known results for the single-component equation 

h t + {h- l h x ) x + (h 3 h xxx ) x = 0, 

which touches down to zero in finite time. This equation has been studied by Burelbach et 
al. [ID], with the initial condition 



1 + 0.1 sin ( x 



They observed finite-time rupture and estimated the the rupture time as t^ = 4.164, based 
on a numerical study with At = 10 -5 and N = 40. With these grid parameters, the 
numerical scheme (A4)-(A5) gives a rupture time t^ = 4.145. A possible source for the 
small discrepancy of estimates is the fact that we have kept At = 10~ 5 for the duration of 
the simulation; Burelbach et al. refine it as rupture approaches, until At = 10~ 5 . Refining 
the grid (N = 100) gives a reduced rupture time tn = 4.093 (see Fig. P), a reduction that is 



consistent with Fig. 3 of Burelbach et al. In conclusion, this test of our scheme validates its 



applicability to the two-component equations (Al) 
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